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It has long been agreed by academics that the inversion method is the method of choice for generat- 
ing random variates, given the availabihty of the quantile function. However for several probability 
distributions arising in practice a satisfactory method of approximating these functions is not avail- 
pi^ I able. The main focus of this paper will be to develop Taylor and asymptotic series expansions 
for the quantile functions belonging to the following probability distributions; Variance Gamma, 
C . Generalized Inverse Gaussian, Hyperbolic and a-Stable. As a secondary matter, based on these 
^+H I analytic expressions we briefly investigate the problem of approximating the quantile function. 
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f^ ■ 1 Introduction 

o' 

^ [ Analytic expressions for quantile functions, have long been sought after. The importance 

of these functions comes from their widespread use in applications of statistics, probability 
theory, finance and econometrics. Therefore much effort has been devoted into their study, 

r> ■ in particular since closed form expressions for the quantile function of most distributions are 

C^ , not known, several approximations appear in the literature. These approximations gener- 

ally fall into one of four categories, series expansions, functional approximations, numerical 
algorithms or closed form expressions written in terms of a quantile function of another 
distribution. The focus of this report is on the former two categories. 

For a given distribution function F denote by Q the associated quantile function. We 
restrict our attention to the class of distributions for which F is strictly increasing and 
absolutely continuous. In this case we have 

Q (u) := F-' (u) (1.1) 

where F~^ is the compositional inverse of F. Suppose the corresponding density function 
/ (x) is known. Differentiating (11. ip we obtain the first order quantile equation, 
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dQju) ^ 1 ^ 1 

du F'{F-Hu)) f{Q{u)y ^-^ 

This is an autonomous equation in which the nonhnear terms are introduced through the 
density function /. Many distributions of interest have comphcated densities, for example 
the densities of the generahzed hyperbohc distributions are written in terms of higher tran- 
scendental functions. Thus the solutions to (11. 2p are often difficult to find, and hence this 
route has been relatively unexplored in the literature. However provided the reciprocal of the 
density 1// is an analytic function at the initial condition Xq = Q (uq) we may employ some 
of the oldest methods of numerical integration of initial value problems, namely the method 
of undetermined coefficients and the method of successive differentiation, used to determine 
the Taylor series expansions of Q. Since the equations in question are non-linear, finding 
their solutions requires some special series manipulation techniques found in for example 
[19j . Shaw and Steinbrecher [42j considered power series solutions to a nonlinear second or- 
der differential equation focusing on particular distributions belonging to the Pearson family. 
We will extend their work here by examining some non-Pearson distributions, in particular 
we will look at the cases where / is the density function of the hyperbolic, variance gamma, 
generalized inverse Gaussian and a-stable distributions. 

The central issue in developing algorithms to approximate the quantile function Q is that 
near the singular points the condition number k,q (u) defined by. 



i^q{u) 



uQ' {u) 



Q{u) 

is large, resulting in an ill conditioned problem (i.e. near the singular points a small change 
in u will result in large deviations in Q). A remedy to this problem is to introduce a new 
variable z to reduce kq. As an example consider the quantile function of the standard 
normal distribution $^^(u) := A/2erf~^ {2u — 1). The associated condition number k^-i (m) 
grows without bound as -u — )■ 1. Hence our problem is ill conditioned near the singular point 
u = 1, and one should therefore seek a change of variable to reduce k$-i (m), for example 
by introducing the variable z (u) := — ln(l — m). One may then consider the equivalent 
problem of approximating the transformation A{z) := $~^(1 — e^^). This approach was first 
taken in [10], albeit from a slightly different perspective. Of course in the normal case the 
transformed function A can be written down explicitly, since the form of Q was know. For 
situations in which this is not possible we formulate a first order differential equation, which 
we will solve in the subsequent sections for A. 

To motivate the idea suppose the base distribution used in the inversion method to gen- 
erate random variates is the standard uniform f/(0,l). Once a random variate from this 
distribution has been generated we may generate a random variate X from the target distri- 
bution with c.d.f. Ft by setting X = Qt (U), where Qt is the associated quantile function 
of the target distribution and f/ ~ f/ (0, 1). However as noted by Shaw and Brickman in [40j 
the base distribution need not be restricted to the standard uniform. Suppose we wish to 
generate random variates from a target distribution with distribution function F^ and asso- 
ciated quantile function Qt- Assume further that rather than generating uniform random 
variates we can generate variates from a base distribution with distribution function Fb and 
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associated quantile function Qb- Now if we could find a function A which associates the u*^ 
quantile from the base distribution with the m*'' quantile from the target distribution, 

Qt{u) = A{Qb{u)), (1.3) 

we could generate random variates from the target distribution simply by setting X = 
^{Qb{U)). Note however the function argument Z = Qb{U) is itself a random variate 
from the base distribution. Thus we now have a recipe to generate random variates from a 
specified target distribution given random variates from a non-uniform base distribution. 

There is no restriction on the choice of base distribution, however for the purpose of 
approximating Qt a wise choice of the variable z := Qb (u) would be one that reduces the 
condition number, ka {z (u)) < Kq (u). We will now derive an ordinary differential equation 
describing the function A arising from the change of variable. Starting with the first order 
quantile equation for the target quantile function Qt, 

(JQt ju) _ 1 

du fAQrHY ^-^ 

we make a change of variable z = Qb {u). A simple application of the chain rule for differ- 
entiation to (11. 3p and using the fact that Qb also satisfies the first order quantile equation 
gives, 

dQriu) dAdz dAdQB dA 1 dA 1 



du dz du dz du dz Jb^Qb (u)) dz fB{,z) 

Substituting this identity into (ll.4p we obtain a first order nonlinear differential equation, 

dA ^ fBJz) 

dz MA{z)y ^ ■ ' 

which we call the first order recycling equation. A second order version of this equation was 
treated in [40j. 

Note that the idea of expressing target quantiles in terms of base quantiles is not a new 
one, indeed this is the idea behind the generalized Cornish Fisher expansion, which originally 
introduced in |11) and [1^ was generalized by Hill and Davis in [23]. The Cornish Fisher 
expansion has many drawbacks however, we refer the interested reader to [2^. Another 
interesting idea along these lines was introduces by Takemura [53] . Here the the Fourier 
series expansion of the target quantile Qt is developed with respect to an orthonormal 
basis of the form {ipi o Qb}^o, where {ipi}°lQ is itself an orthonormal basis for a set of 
square integrable functions, we refer the reader to the original paper for details. Unlike 
the Cornish Fisher expansion which is asymptotic in nature Takemura's approach yields a 
convergent series in the L^ norm. Note however the computation of the Fourier coefficients 
usually requires numerical quadrature and that for approximation purposes the L°° norm is 
preferred. 

Standard numerical techniques such as root finding and interpolation for approximating 
the associated quantile functions often fail, particularly in the tail regions of the distribution, 
see for example [14j. Thus in addition to developing convergent power series solutions to 
(II. 2p and (II. 3p we also develop asymptotic expansions of Q at the singular points u = 
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and u = 1. In general it is not easy to discover these asymptotic behaviors, some trial and 
error and intelligent fiddling is required, [51 Ch. 4]. Our approach is to consider Q as being 
implicitly defined by the equation, 

F(Q) = u. 

Asymptotic iteration methods, see for example |i8j, Ch. 2], may then be used to obtain 
the leading terms in the expansion. For numerical purposes the resulting expansions are 
divergent but may be summed by employing an appropriate summation technique such as 
Shank's e^ or Levin's u transforms. The reader is warned that the approach taken throughout 
the paper is an applied one, the reader will certainly notice a lack of rigor at certain points, in 
particular in deriving these asymptotic expansions we will execute some "sleazy" maneuvers. 
In this report we will focus on the variance gamma, generalized inverse Gaussian, hyper- 
bolic, normal inverse Gaussian and a-stable distributions. For each of these distributions we 
will find power series expansions of Q and A as well as the asymptotic behavior of Q near its 
singularities. Despite the unsightly appearance of the formulae given in the present paper 
for the coefficients appearing in these expansions, they are simple to program if not tedious. 
It is difficult to study analytically the convergence properties of these series, hence we will 
proceed in an empirical manner. In the last section of this report we experiment with various 
numerical algorithms to approximate Q, observing that a scheme based on Chebyshev-Pade 
approximants and optimally truncated asymptotic expansions seems most useful. 

2 Hyperbolic Distribution 

Amongst others it was observed by [16] and [6] the hyperbolic distribution is superior to 
the normal distribution in modeling log returns of share prices. Numerical inversion of 
the hyperbolic distribution function was considered in [25]. There numerical methods were 
considered to solve (11. 2p and only the leading order behavior of the left and right tails was 
given. Here we will provide an analytic solution to (11. 2p and the full asymptotic behavior of 
the tails. The density of the of the hyperbolic distribution Hyp [a, /3, 6, fi) is given by, 

/ (x; a, (3, 6, /i) = ^ ^ J ,^ -a^s^H^-,)--^^i^-,)^ (2.1) 

2adKi (()7) 

where a > 0, \(3\ < a, 6 > and ^ ^M. are shape, skewness, scale and location parameters 
respectively and for notational convenience we have set 7 = \/a^ — /3^. By defining ai := 5a 
and Pi := 5/3, we obtain an alternative parametrization in which ai and Pi are now location 
and scale invariant parameters. Hence without loss of generality we may set yU = and 6 = 1, 
since Q {u; ai, Pi, fi,6) = fi + 6Q {u; ai. Pi, 0, 1). The first order quantile equation (ll.2p then 
reads, 

^ = Noe'^rV^^-^iQ (2.2) 

du 

where A'"o = 2aiKi (71) /71 and 71 = ■\/a;f — Pf. To form an initial value problem (IVP) let 
Uq e (0, 1) and impose the initial condition Q (mq) = Xq. For all practical purposes xq is 
usually determined by solving the equation F (xq) — Uq = using a root finding procedure. 
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By applying the method of undetermined coefficients we find Q admits the Taylor series 
expansion, 

oo 
n=0 

where the coefficients g„ are defined recursively as follows, 



Xo, n = 

'"^^^&.-i, n>l 



\Yl=iK(^io.k - l^iqk)hn-k, n>l 



and, 



a/1 + Xq n = 

;^ (^gngo + YJiZli^ + i)(gfc+ign-fc-i - 0^+10^-^-1)) n>i 

To develop the asymptotic behavior at the singular points -u = and u = 1 note that the 
hyperbolic distribution function F satisfies the relationship, 

F (x; ai, /3i, 1, 0) = 1 - F (-x; ai, -/3i, 1, 0) , 

which implies, 

Q {u; ai, /3i, 1, 0) = -Q (l - u; Oi, -/3i, 1, 0) . 

Hence without loss of generality we need only look for an asymptotic expansion of Q as 
M — > 1. Now under the assumption 1 ^ Q as m — ;■ 1, we obtain from fl2.2p . 

du 
Solving this asymptotic relationship along with the condition Q {1) = oo yields the leading 
order behavior of Q, 

Q[u)^ , asM-)-l. 

(ai - I3i) 

From which we note that Q has a logarithmic singularity at m = 1. To develop further 
terms in the expansion, write x := Q (u) from which it follows u = F{x) and note that Q is 
implicitly defined by the equation. 



1-F{x)= f{t)dt 

J X 

Rearranging and introducing the variable v we obtain. 
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/■oo 

v:={l-u)No= / e-"^^^+*^+^^*c^t 



The idea is to expand the integrand appearing in the right hand side and integrate term- wise, 
this process can be carried out symbohcally, giving the first few terms in the expansion, 

\ai-/3i 2{ai-^i)x 8{ai - Pi^x^ \x^ 

Taking logs and inverting the resulting series allows us to write down the first few terms in 
the asymptotic expansion of x, 

«i Qfi n f ^ 

from which we can conjecture the form the asymptotic expansion of Q as, 

oo 

Q{u)r^y + y^^, asu^l, (2.3) 



where, 

ln(iVo(ai-/3i)(l-n)) 
(ai - f3i) 

Substituting fl2.3p into the first order quantile equation fl2.2p allows us to derive a recurrence 
relationship for the coefficients g„, 

'0 n = 



ai 



where. 



and. 



2{ai-/3i) '^-^ 



n = 

i + gi n = 1 , 
qn + dn n>2 

_ 1 n = 

n SLl ^bkCn-k n > 1 ' 
dn = 7 / (k + l){qk+iqn-k-2 — 0,k+lO'n-k-2] 

n — I ^-^ 

fc=0 
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Note that (12. 3p is a divergent series, however of the sumination methods we tested we 
found Levin's u transform and Pade approximants particular useful for summing (12. 3p . In 
both cases analytic continuation was observed. Later we will briefly look at algorithms 
for constructing rational approximants of Q valid on the domain [10~^*^, 1 — 10~^°], but if 
necessary one may utilize (12. 3p to obtain approximations on a wider region. 

As mentioned in the introduction the problem of approximating the quantile function 
near its singular points is an ill conditioned one. To this end it will be useful to introduce 
a change of variable which reduces the condition number hq. Motivated by the asymptotic 
behavior of Q near its singular points, in particular its leading order behavior we introduce 
the base distribution defined by the density. 

Here Xm '■= /3i/7i is the mode of the hyperbolic distribution, p_ := e~^°'^~^^^^^"^Pm, p+ '■ = 
g(ai-/3i)xm ^i — p^) and Pm = Fx{xm)- The associated distribution and quantile functions 
can be written down as, 
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Fb ix) := < , a \ 

|l-p+e-("i-^i)^, X>Xr, 

and, 

Qb{u) 

respectively. Substituting this choice of fs into the recycling equation (II. 5p results in a left 
and right problem, 

dz 
and 

dA 



^ln(i^), u>p. 



dz 
respectively, along with the suitably imposed initial conditions. For the left problem, we 
choose Mo G {0,Pm] and impose the initial condition Xq = A{zq) = Qt{uq), where Zq := 
Qb {uq)- Similarly for the right problem choose mq G [p_, 1). Again through the method of 
undetermined coefficients we obtain the Taylor series expansion of A, 



^(^) = "^a^iz- zo) 



n=0 



where the coefficients are defined recursively by. 
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and. 



Xq n = 

^e {ai + Pi) dn-i n>l 



'l + xt. 



n = 



^ (nantto + I]fc=o(^ + 1) {ak+ian-k-i - bk+ibn-k-i)) n > 1 
"i {h + {/)) + /3i (1 - fli) n = 1 



d„ 



n Y2=l kCkdn~k n>l 

For the solution to the left and right problem make the following replacements, 



Tab. 1: Coefficient Parameters 





Left 


Right 


e 


P~ 


P+ 





+1 


-1 


P 


(«i + A) 


- (ai - /3i) 



3 Variance Gamma 

The variance gamma distribution was introduced in the finance literature by Madan and 
Seneta in [30]. To our knowledge there has been very little written on the approximation of 
the variance gamma quantile. The density is given by, 



/vG(a;;A,a,/3,/i) 
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2A 



|A-i 



vA-1/2 



{2ar-'^'y^T{X) 



X — /i| 2 K^_i (a |a; — /i|) e 



/3(x'-m) 



(3.1^ 



where A > 0, a > 0, |/3| < a, 7 = y/a — (5 and ^ &^. Setting the location parameter /i to 
zero and substituting the density (13. ip into the quantile equation (11. 2p gives. 



du 



(3.2) 



where the function g is defined as. 



9 {y) ■■= No 



e-/3?;|^|f-A 
KA-i("|y|)' 
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and A^o = (2a) y^^ (A) /7^^. Our strategy to solve (13. 2 p will be to apply the method 

of successive differentiation to obtain the Taylor series representation of Q, 



Qyg (u) = ^qniu- Mo)" 



(3.3) 



n=0 



where go is determined by the imposed initial condition, and the remaining coefficients are 
given by, 



d 



n-l 



Qn 



du 



n-l 



g{Q{u)), n>l 



Thus the problem reduces to finding the higher order derivatives of the composition g o Q^ 
which can be obtained recursively through by applying Faa di Bruno's formula as follows. 
Note first that g can be written as the product of three functions a [y) := e~^^, b {y) := |y|2^'^ 



-1 



and c (y) := K^_i{a\y\) . Thus an apphcation of the general Leibniz rule yields. 



n k 



^^"^ (y) = A^o E E it) {^) ^'^'"' (^) ^''' (^) ^''"' (^) • 

fc=o i=o V / V-^/ 
The higher order derivatives of the functions a and h appearing in (13. 4p are given by. 



(3.4) 



«(") (y) = (-/?)" e- 



-Py 



feW {y) 



y\-^-^ 



n-l . 

n i 



1=0 



A — 2 X 



-1 i/>0 
-1 y<0 



To find the n*'* derivative of c note that c can be written as the composition c (y) = Ci 0C2 (y), 
where Ci (y) := Xjy and 02(1/) := K y^_\{pi\y\) . Hence we may obtain c^") (y) through an 
application of Faa di Bruno's formula. 



.(") 



M^E 



n\ 



milm2l ■ ■ -nir, 



tC 



(miH \-m„) 



.C2{y))l[ 

j=0 



SJ) 



cT (y) 



j! 



where the summation is taken over all solutions {1711,7712, . . . , 7nn) € Z>q to the Diophantine 
equation. 



7711 + 2r7i2 + 3?7l3 ■ ■ ■ + 7l7nn = 71. 

The formulae for the higher order derivatives of Ci and C2 are given by, 
and 



(3.5) 



Jn) 



or iv) 



2 



E 

A;=0 



n\ ^ iK;,_i_(2fc_„)(ay), 



2/>0 
-ay) y < 



3 Variance Gamma 10 

Where we have used the identity [33l eq. 10.29.5], 

Kl") (z) = (-1) " fl (l) K-(2/.-n) (z) , (3.6) 



^ ' k=0 ^ ' 

which can easily be proved by induction. Formula (13.61) shows that, to find the higher order 
derivatives with respect to z of the modified Bessel function of the second kind K^"^ (z) we 
need only a routine to compute Kj, [z). Now given the scheme (13. 4p to compute g^^' {x) 
and the coefficient go defined by the initial condition we may then compute g„ for n > 1 
recursively, by another application of Faa di Bruno's formula. It is important to note that 
to generate the coefficients g„ appearing in (13. 3 p does not require the use of any symbolic 
computation. 

Next we will focus on deriving an asymptotic expansion for Q (u). Similar to the hyper- 
bolic quantile the following equality holds, 

Q [u] A, a, /3, 0) = —Q (1 — u; A, a, — /3, 0) , 

so again without loss of generality we need only seek an asymptotic expansion of Q as m — )■ 1. 
Our strategy here will be to, 

1. derive an asymptotic expansion for the density / as a; — > cxd, 

2. integrate term-wise to obtain an asymptotic expansion for the distribution function F 

as s — 7- oo, 

3. and finally invert this expansion to obtain an asymptotic expansion for the quantile 
function Q as m — )■ 1. 

For the first step we will make use of the asymptotic relationship [33| §10.40], 

1 oo / X 



A;=0 



where, 

k 

k\8~ 



«^(^)^=T4^n(4-'-(2j-in- 



From which it follows 



(2«)'r(A) ^ «^ 

Working under the assumption that term-wise integration is a legal operation we obtain an 
asymptotic expansion for the distribution function. 
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2A °^ /"oo 



l-F(x) ~ 7^^,y2<y~'~'a,{\-2~') I e-^^-f^y^t'-'^-'dt 



2A ^^ 

= ^i|77TyE""'''«'^(^-2"')r(^-fc'^(«-/3)), (3.7) 

where F (a, z) is the upper incomplete gamma function, which for large z satisfies the asymp- 
totic relationship [33| §8.11], 

^r(a-j) 
Substituting into (13. 7p and assuming the terms of the series may be rearranged we obtain, 

/o \^^ 2A ''-^ 

l-F(x)^ ^ "^1.^^ x^-^e""^"-^) V&fcX~' asa;^oo, (3.9) 

where, 

6, = ^^(a - /3)-(^+i)a-('=-^) ( JJlA - fc + ^) ) a,., (a - ^) • 
j=0 \j=o / ^ ^ 

The expression 03. 9p describes the asymptotic behavior of the variance gamma distribution 
function as x — )■ oo. Let u = F (x) in (13. 9p . our goal then is to invert this relationship to 
obtain an asymptotic expansion of the quantile function Q (m) as m — )■ 1. Introducing the 
variable, 

\A 



and rearranging (13. 9 P we obtain 



^^ = yA 1^-^^ 



w~a;^"^e~^("~^)D(ij as x ^ oo. (3.10) 

where D is the formal power series defined by D (z) = XlfeLo ^kz''- Taking logs and introduc- 
ing the variable y := — In w/ {a — /3), we may write (I3.10p as, 

A-1, 1 , fl\ 

x^y + ^logx + -logD - (3.11) 

a — [3 a — [3 \x J 

We wish to write x in terms y; this task might at first may appear difficult to achieve, 
but as it happens we are in luck, similar expressions occur frequently in analytic number 
theory, and some useful methods have been developed to invert these kinds of relationships. 
A drawback of these methods is that they require symbolic computation. The most basic 
method known to us, one can apply to invert (13. lip is the method of asymptotic iteration 
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[8]. This method however is extremely slow. A much more efficient approach is to use the 
method of Salvy[T3]. There it was noted that the form of the asymptotic inverse is given by, 

n=0 ^ 

where Pq (^) is a polynomial of degree 1 and P„ (^) are polynomials of degree n for n > 1. 
Following a similar analysis of that in [13] it can by shown that Pn may be determined up 
to some unknown constant terms cq, . . . , c„ by the recurrence relationship, 






^"^^^ ^M(^n-i(O-i^n-i(0))-|^(n-l)/„«P„_i(t)rft + c. n>0- ^^-^^^ 



The unknown terms cq, . . . , c„ may be computed through the following iteration scheme, 

• starting with uq (t) = log {bo), compute, 

Uk (t) = ; 'i\ log (1 + tuk-i (t)) + -^ log D (-—-^ ^ ] , k = l,...,n + l, 

• extract the constants, 

Ck = [t''] Un+i (t) , k = 0,...,n 

where [t*^] m„_i-i (t) is used to denote the coefficient of the t'^ term in Un+i (t). 

An implementation of this scheme in Mathematica code has been included in appendix 
(TODO: Include code). Through this process the first few terms of the asymptotic expansion 
of Q may be generated as follows, 

^, . , (A - 1) logy -log (a -/3) 
Q{u) ~ y + ^ 

a — p 

(A - 1) (2a + aX- PX + 2a{X - 1) logy - 2a log [a - /?)) 
+ 2ya{a - Pf +""'' ^^ ^ "^ ^^ 

Next we focus on introducing a change of variable to reduce the condition number hq 
near the tails and finding the Taylor series representation of the corresponding function A. 
As in the hyperbolic case motivated by the asymptotic behavior of the quantile function near 
its singularities, in particular by its leading order behavior we introduce the base distribution 
defined by the density, 

_ J]9_ (a + /3) e('^+^)^ X < 
^"^ ^""^ •" \p+ {a - 13) e-(-«- X > ' 

where p_ = F (0) and p+ = 1 — p_. The associated distribution and quantile functions can 
be written down as, 
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P-e 



(ai+/3i)x ^ <; g 



and 



Fb ix) :-— . , a\ 



-^ In H^ , u<p. 



respectively. Substituting this choice of fs into the recychng equation fll.Sp results in a left 
and right problem, 



dA , , ^ , 

— =gL{z), z<Qb{P- 



and 



dA 

-r = 9Riz), z>Qb{p^), 
dz 

respectively, along with the suitably imposed initial conditions. For the left problem, we 

choose Uq G (0,p_] and impose the initial condition Xq = ^ (^o) = Qt{uo), where Zq : = 

Qb {uq). Similarly for the right problem choose Uq G [p-, 1). The functions Ql and Qr 

appearing on right hand side of these differential equations are defined as, 

g^{z):=p.{a + P)e^^^^>g{A{z)), 
and, 

gn{z):=p^{a-P)e-^^-f'>g{A{z)). 
Suppose that the series solution of either problem is given by, 

oo 

A{z) = J2(^n{z~Zor. 
n=0 

Here the first coefficient Oq is determined by the initial condition imposed at Zq and the 
remaining coefficients are given by, 

an = Ql'^^ (^o) , n>l, 
for the left problem and, 

an = 9k~^^ (^o) , ri>l, 

for the right problem. Both sets of coefficients may easily be computed from an application 
of Liebniz's rule, 



^[") (z) = p_e(°+«^ 



j:[f)ic. + firHgoAt-'Hz), 
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Starting with oq, the higher order derivatives of the composition g o A appearing in these 
formulae, are computed recursively in precisely the same way we computed the higher order 
derivatives oi g o Q above. 

4 Generalized Inverse Gaussian 

The generalized inverse Gaussian (GIG) is a three parameter distribution, special cases of 
which are the inverse Gaussian, positive hyperbolic and Levy distributions to name a few. It 
arises naturally in the context of first passage times of a diffusion process. The probability 
density function of a GIG random variable is given by, 

/gig(x;A,x,^) = (^y ^ A-i^^Kx^-^+V^^), ^>o, (4.1) 

where AgM, x>0,'?/^>0 and K^, (z) is the modified Bessel function of the third kind 
with index v. The GIG distribution is also used in the construction of an important family 
of distributions called the generalized hyperbolic distributions; more specifically a normal 
mean mixture distribution where the mixing distribution is the GIG distribution results in 
a generalized hyperbolic distribution. Consequently if one can generate GIG random vari- 
ates then a simple transformation may be applied to generate variates from the generalized 
hyperbolic distribution [48]. 

We will use an alternative parametrization to the standard one above, let i] = a/x/V and 
i^ = Vxi^j the density then reads. 

In this new parametrization u and A are scale invariant and 77 is a scale parameter, so in 
the following without loss of generality we may set r; = 1. The first order quantile equation 
(II. 2p now reads, 

^ = 2KA(a;)eH^+Q)Qi-\ 
du 

Let uq G (0, 1) and impose the initial condition Q (uq) = xq. For the case A 7^ 1, the GIG 
quantile Q admits the Taylor series expansion, 

00 
Q{u) = J2(ln{u-uor, (4.3) 

n=0 

where the coefficients g„ are defined recursively as follows, 

Xq n = 

' ^Kxiuj) J2'^^Q hcn-i-i n > 1 ' 
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^ n = 

10 

1 \~^n 



- Er=l ?i«n-i "^ > 1 



<?0 



, 2n 

and 






gi-^ n = 

1 Y^Ti /(2— A)« 



«o 



i:ELi ^^-lU.c_ n>l 



For the special case A = 1, the coefficients are somewhat siniphfied, with a^ and 6„ as above 
the coefficients appearing in ( 14. 3 p become, 



Qn , 2 



Xq n = 

lKi{uo)hn-i n>l 



Next we will focus on developing the asymptotic behavior of Q as m — )■ 1. We proceed in 
an analogous fashion to the variance gamma case, and find that remarkably the form of 
asymptotic expansion of Qgig is very similar to that of QvG as m — )■ 1. From the definition 
of the distribution function we have, 



oo 



Expanding the e""^/^* term and integrating term-wise we obtain, 

'^^ ^ fc=0 ^ ^ 

where F (a, z) is the upper incomplete gamma function. As in the variance gamma case 
substituting (13.81) and rearranging the terms provides us with a more convenient form of the 
asymptotic expansion for the generalized inverse Gaussian distribution function, 

^ oo 

l-F(x)^ x^-^e-f^bkX-'', asx^oo, (4.4) 



2Kx(uj, 



where 



''=E^(i) n(^-*+' 

Now let M = F (x) and introduce the variable v = 2Kx{u)){l — u), then we can rewrite (14.41) 



as, 



vr^x^-^e--D(-), 4.5 

X 
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where D is the formal power series defined by D {z) = 'Yl^=o^kZ^ ■ To invert the asymptotic 
relationship (14. 5 p we start by taking logs and introducing the variable y := —(2/uj) Inf. We 
may now write (14. 5p as, 

X r-u y -\ log X -\ — log D { — ] . (4.6) 

U UJ \x J 

One may now apply the method of asymptotic iteration [8j to invert (14. 6p . As mentioned 
earlier this method however is extremely slow and a much more efficient approach is to use 
the method of SalvyfTS]. Again the form of the asymptotic inverse is given by, 



x = Q{u) ^y + 2_^ 



yn 
n=0 ^ 



where Pq (^) is a polynomial of degree 1 and P„ (^) are polynomials of degree n for n > 1 
and following a similar analysis of that in [13j it can by shown that P„ may be determined 
up to some unknown constant terms cq, . . . , c„ by the recurrence relationship. 



-^^i + Co n = 

^(Pn-l (0 - Pn-1 (0)) - ^(n - 1) /,^P._i (t) dt + Cr. n>0 



The unknown terms Cq, . . . , c„ may be computed through the following iteration scheme, 

• starting with uq (t) = log (bo), compute, 

Uk{t)= ~ M og(l + fafc_i(t)) + -logD(— — ^ --) , k = l,...,n + l, 

• extract the constants, 

Cfc = [t^] Un+l (t) , /c = 0, . . . , n 

where \t^'\ Un+i (t) is used to denote the coefficient of the t'' term in Un+i (t). 

Through this process the first few terms of the asymptotic expansion of Q may be generated 
as follows, 

2((-l + X)lny + ln(^)) 
Q{u) ~ y+^ ^^ ^-^ 

4A - ^2 _ 4 ^ 4(_X ^ A)2 Ini/ + 4(-l + A) In (^) 

H ^^ H , as u -)■ 1. 

u^y 

To observe the asymptotic behavior of Q (u) as m — )■ 0, we utilize the following identity, 

Q (w; A, 1, uj) = -— \ ^ . , (4.8) 

Q{1- u] -X,l,uj) 



which can be easily proven as follows; by definition of the density (14.21) we have. 
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/ (x; A, 1, uj) = —f i -; -A, 1, w ) . 



x^' \ a; 



Integrating both sides then yields, 



F(x;X,l,u) = 1-F ( -;-A, l,w 

\x 

from which (14. 8 P follows. 

Next we consider solving the recycling ODE, but first we must choose a base distribution. 
Again motivated by the asymptotic behavior of Q as m — )■ and -u — )■ 1, in particular 
the leading order behaviors we suggest the following base distribution characterized by the 
density function, 

fB i^) = { :f _-. . • (4-9) 



UJ „-^x 



PRI^ 2 



, X ^ *^?7l 



Here Xm serves as a cutoff point between two suitably weighted density functions, in partic- 
ular Xm is the mode of the GIG distribution defined by. 



UJ 

Note that the density function fi [x] := (uj/2x'^)e^'^^'^^ belongs to the scaled inverse x^ 
distribution with 2 degrees of freedom and scale parameter u/2, and that the density function 
fji (x) := (1^/2) e""^^/^ is the density of an exponential distribution with rate parameter uj/2. 
The normalizing constants pl and p^ are defined by, 

PL = e^Pm, 

and, 

PR = et^'" (1 - p„) , 

where Pm = Fqiq (xm) ■ The associated distribution and quantile functions can be written 
down as, 

p f X J Plc"^, X <x„^ 

Fb {x) = < ^^, . 

l-(l-pR)e 2^^ ^"^ x> x,n 



and. 



U<Vn 



Qb{u) = { ^^^("/^,-) ^ ,. "-"'", (4.10) 



respectively. Substituting this choice of Jb into the recycling equation (II .Sp then leads to a 
left and right problem given by. 
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dA 
and, 



^^ 2plKa {uj) ^e-^e^-(^+^)Ai-\ (4.11^ 



^=PRu;K,MeHi+^-)A^-\ (4.12) 

respectively, along with the suitably imposed initial conditions. For the left problem, we 
choose Uq G (0,Pm] and impose the initial condition Xq = A{zq) = Qt{uo), where Zq := 
Qb {uq). Similarly for the right problem we choose Uq G \pm, !)• Treating the left problem 
first, we find, 

oo 

A{z) = Y,(^n{z-zor, (4.13) 

n=0 

where the coefficients are computed recursively through the identity, 

Xq n = 



ar. 



.^K, (c) (e::J E;:o"'(^ + ^)h^id,e^^k-,-: 



n > 1 ' 



where, 









and 



^ n = 

ao 
gtj(co+ao)/2 n = 



ag"^ ra = 

'"='iEL(^-l)a.e„-. «>1 



The coefficients appearing in the series solution f l4.13p to the right problem 14.121 are given 
by. 



Xo n = 

.^K,(u;)(n:;4e„-fc-i) n>l' 

where c„ and e„ are defined as in the solution to the left problem above and, 

■ 2(01 + 01-1) n = l 
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and, 

(^bo~ujzo/2 ^ ^ Q 

\^Yl'i=i'^bidn-i n > 1 
5 a-Stable 

Under the appropriate conditions Lagrange's inversion formula is capable of providing us with 
a series representation of functional inverses. Yet it seems to be ignored in the literature 
when one wants to find an approximation to the quantile function, which itself is at least in 
the continuous case, defined as the functional inverse of the c.d.f. Based on this observation 
we provide a convergent series representation for the quantile function of the asymmetric a- 
stable distribution. Note however the method is much more general than this; all it requires 
is a power series representation of the c.d.f. 

Suppose the c.d.f. Fx and quantile function Qx of a random variable X have the Taylor 
series representations, 

oo „ oo 

Fx (x) = y^^{x - xo)" , and Qx (u) = y^^{u- uq)" , 

71=1 71=1 

respectively. The relationship between Fx and Qx is given by Fx {Qx (u)) = u, for all 
u G (0, 1). The goal is to solve this expression, that is we would like to write the Taylor 
coefficients g„ in terms of /„. One such expression to achieve this is provided by Lagrange's 
Inversion formula written in terms of Bell Polynomials, see [21 § 13.3], 

jfi ^= ^ 

[~7^ Z]fc=l (lk^n,k (/l, • • • , /n-fc+l) n> 2 

where the coefficients Mn^k (/i, • • • , fn-k+i) are the Bell polynomials. They are defined by a 
rather detailed expression, 

^n,k[jl, ■ ■ ■ , Jn-k+l) '■= 2_^ yTE r , / -U^j] A J2 " ' Jn-k+1 ^ (^'2) 



vi+V2+---=k 
Vi+2v2+3v3---=n 

where the summation is taken over all solutions (fi,f2, . . . ,f„) G Z>q to the Diophantine 
equation, 

vi + 2f 2 + 3t;3 h nvn = n (5.3) 

with the added constraint, the sum of the solutions is equal to k, i-e.^"^Qfj = k. For 
example the solutions (f i, f2, fs, v^ for n = 4 are given by 

{(4,0,0,0), (2,1,0,0), (0,2,0,0), (0,0,0,1), (1,0,1,0)}. 

and if A; = 2 this picks out the solutions (0,2,0,0) and (1,0,1,0). Note that solutions to 
(15. 3 p correspond exactly to the integer partitions of n. An integer partition of a number n 
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is an unordered sequence of positive integers who's sum is equal to n. The added constraint 
imphes we should look for partitions in which the number of non zero summands is equal to 
k. For example the integer partitions of n = 4 are given by (1, 1, 1, 1), (1, 1, 2), (2, 2), (4), 
and (1,3), and in the case k = 2 this singles out the partitions (2,2) and (1,3). Thus to 
summarize the sum in (15.21) is taken over all integer partitions of n in which the number of 
summands is given by k. 

The above recursion (15.1 p can be solved (see [?, § 13.3]) leading to a more computationally 
efficient expression for the coefficients, 




n = 1 



'^> (n-1)! ^"-l.fc (^2/i' 3/i' • • • ' (n-fc+l)/i j n>2 

The a-stable distribution, denoted S^ (/3, /x, cr) is commonly characterized through its 
characteristic function. However there are many parametrizations which have lead to much 
confusion. Thus from the outset we state explicitly the three parametrizations we will 
work with in this report and provide the relationships between them. We will call these 
parametrizationsPo, Pi and P2 respectively. In the following a subscript under a parameter 
denotes the parametrization being used. We will primarily work with Zoltarev's type (B) 
parametrization, see [50^ p. 12] denoted by P2. In this parametrization the characteristic 
function takes the form, 

(t) = exp {a, {ttfi, - \tr e-^(-/2)fti^(-)«g°(*)) } , (5.5) 

where a G (0, 2] is the tail index, /i2 G K is a location parameter, o"2 > is a scale parameter, 
P G [—1, 1] is an asymmetry parameter and K{a) := a — 1 + sgn(l — a). Pi is the classic 
parametrization, and is probably the most common due to the simplicity of the characteristic 
function given by, see [39| p. 5], 

(t) = exp I -< |t|" (1 - i/Sitan (—) sgn (t)) + ifiitl . (5.6) 

A connection between the parametrizations P2 and Pi can be derived by taking logarithms 
and equating first the real parts of (15. 5p and (15. 6p . followed by the coefficients of t and |t|" 
in the imaginary parts, leading to the set of relations. 



Hi = H2cr2 

(Ti = ( cos ( -ttK (a) (32 j 02 

/7ra\ /I 
/3i = cot (^— j tan ( -'kK {a) f32 

Despite the fact no closed form expression for the c.d.f. Fa {x; /3, fi, a) of the stable distribu- 
tion in the general case is known, it can be expressed in terms of an infinite series expansion. 
As is usual for location-scale families of distribution, without loss of generality we may set 
the location fi2 and scale (T2 parameters to and 1 respectively. In addition it is sufficient 
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to consider expansions of the c.d.f. F^ (x; /3, 0, 1) for values of x > only since the following 
equality holds, 

F„ (x; /3, 0, 1) = 1 - F„ {-X- -/3, 0, 1) . (5.7) 

Theorem 1. Let X ~ ^^ (x; /32, 0, 1) he a standard a-stable random variable. Then the 
cumulative distribution function of X admits the following infinite series representations, 



n=0 



F^{x; 132,0,1) = f2^x\ (5. 



where, 

In 

and 




l32K{a) 



n = 
n>l 



F„(a;;/32,0,1) = V^x-"", (5.9) 



n! 

n=0 



~^fl n = 

If 1 < a < 2 then Ii5.8\) is absolutely convergent for all x > and if < a < 1 then Ii5.8\) 
is an asymptotic expansion o/ F^ (x; /32, 0, 1) as x — ;■ 0. In contrast if < a < 1 then the 
series i\5.9\) is absolutely convergent for x > and if 1 < a < 2 then Ii5.9\) is an asymptotic 
expansion of F^ (x; (32, 0,1) as x —)■ oo. 

Proof. The proof proceeds by obtaining first a series expansion of the density of X. This 
is achieved by applying the inverse Fourier transform to (15. 5p . expanding the exponential 
function, and then performing a contour integration, see [2^ § 5.8] for details. The expansion 
in the density can be integrated term by term to obtain an expansion of the c.d.f. D 

Note when 1 < a < 2, the series (15. 8p rapidly converges for values of x near zero, where 
as (15. 9 P converges rapidly for large values of x. The opposite is true when < a < 1. Thus 
we can now apply Lagrange's inversion formula to find a series expansion of the quantile 
function Q^ in the central and tail regions. Note however that since the expansions (15. Sp 
and (15. 9 p are valid only for x > 0, the resulting expansions of Qa are only valid for u > Uq, 
where Uq is the zero quantile location defined by Uq := F^ (0; (32, 0, 1). This does not pose a 
restriction however since it follows from (15. 7p . 

Qa{u; P, 0,1) = ~Qc.{l - u; -P, 0,1) . 

Applying Lagrange's inversion formula to (15. 8 P we see that the quantile function has the 
following infinite series representation valid for u > Uq, 

oo 

Q,(n;/3,0,1) = V^(M-no)", (5.10) 

n=l 
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where the coefficients g„ are given by (15. 4p and 



1/ (32K{a) 

"0 = 7T 1 - 



2 V a 

To find the functional inverse of (15. 9p we make the change of variable y := x~ , and apply 
Lagrange's inversion formula to the power series, 

n=0 

The quantile function is then given by, 

Q„(m;/3,0,1)= [G-'{u)y^, 
where, 

oo ~ 

G"'w = E^ ("-!)"' (5.11) 

n=l 

and the coefficients g„ are given by (15. 4p with /„ replaced by /„. Note when 1 < a < 2, the 
series ( I5.10p rapidly converges for values of u near Mq, where as ( 15. lip converges rapidly for 
values of u close to 1. In this case partial sums of (I5.10p serve as good approximations of 
the quantile function in the central regions where as partial sums of (15. lip can be used to 
approximate the tails. The opposite is true when < a < 1. The first few terms of (I5.10p 
are given by, 

n( R n ^\ 7rcsc(7rp) Tricot (vrp) esc (vrp) T {^^) ^2^c^(( \^\ 

Qa{u; /3, 0, 1) = ^ . , ^ [u - Mo) + —3 [u - Mo) +0 [{u - no) j 

Implementation is a rather straightforward matter, high level languages such as Mathe- 
niatica have a built in implementation of the Bell polynomials and plenty of algorithms exist 
to generate integer partitions in lower level languages such as C++, see for example [35]. 
As long as we have a series representation of the c.d.f. using this approach we could derive 
a series representation for the associated Quantile function. However there is one obvious 
drawback, even though the coefficients g„ can be computed using elementary algebraic op- 
erations, the number of partitions p (n) of an integer n grows exponentially with n, thus for 
large values of n the sum in (15.21) may be computationally expensive due to the large number 
of summands. 

However the problem of reverting a power series is a classical one in mathematics and 
many efficient algorithms have been devised as a result. For example Knuth [271 § 4.7] gives 
several algorithms for power series reversion including a classical algorithm due to Lagrange 
(1768) that requires O [N^) operations to compute the first A^ terms. More recently Brent 
and Kung [7] provide an algorithm which requires only about 150 (A^ log A^) floating point 
operations. Dahlquist et al. ^2j also present a convenient but slightly less efficient algorithm 
based on Toeplitz matrices. 
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Concerning the numerical evaluation of the distribution function F^ of the stable distri- 
bution it has been remarked by various authors [43j that the series expansions given in [T] 
are only useful for approximating F^ for either small or large values of z, due to the slow 
convergence of the series. It is for this reason standard methods such as the Fast Fourier 
transform or numerical quadrature techniques are applied to evaluate Fa- However we found 
in our experimentation that series acceleration techniques such as Fade approximants and 
Levin type transforms [17] could be applied to fl5.8p and fl5.9p . The resulting rational func- 
tions were affective for approximation purposes. The same comments apply for the series 
representation of the density and quantile functions, at least for the set of parameters we 
tested, which include those occurring frequently in financial data. We will discuss numerical 
issues further in the next section. 

6 Numerical Techniques and Examples 

The goal of this section is to discuss the design of an algorithm which accepts a set of 
distribution parameters and constructs at runtime an approximation Qa to the quantile 
function satisfying some prespecified accuracy requirements. In particular we are interested 
in parameters occurring most often in Financial data, but the proposed algorithm works for 
a much wider range of the parameter space. For the distributions we have considered in this 
paper, most published algorithms of this nature are based on interpolating or root finding 
techniques. We break this mold by briefly examining other methods of approximation based 
on certain convergence acceleration techniques. We shall call the time it takes to construct 
an approximation the setup time, and the time it takes to evaluate Qa at a point u the 
execution time of the algorithm. With the analytic expressions made available in this report, 
a wide range of possibilities become available. For instance the expansions may be used in 
conjunction with other numerical techniques: 



• Techniques based on interpolation often fail in the tail regions [15], for this reason 
a good idea would be to supplement the algorithm with the asymptotic expansions 
developed above. 

• Root finding techniques are known to converge slowly. To improve the rate of conver- 
gence one needs to supply a good initial guess of the root. Such as provided by the 
truncated Taylor series provided above, or better yet a Fade approximant. 

Another plausible approach is to construct a numerical integrator based on the Taylor 
method. It has been reported by many authors that when high precision is required this is 
the method of choice, see for example [10], [1] and [25]. The idea here is to discretize the 
domain [a, h] C (0, 1) into a non-uniform grid a = Uq, . . . ,Un = b. To build this grid we must 
determine the step sizes hk = Uk — Uk-i for k = 1, . . .n. In addition at each grid point Uk 
we must determine the order nik of Taylor polynomial so that the required accuracy goals 
are achieved. That is both the stepsize and order are variable. Based on two or more terms 
of the Taylor series certain tests have been devised to compute the hk and mfc, see again the 
references mentioned above and the review article 1201. 
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The result of Taylor's method is a piecewise polynomial approximation to the quantile 
function. However in this report we are more interested in constructing rational function 
approximants of the form, 



Rr, 



Pm (v) _ E^O (^i'"' 



and so we will not discuss Taylor's method further. Some of the best known algorithms 
for approximating quantile functions are based on rational function approximations. Unfor- 
tunately such approximations traditionally are only available for "simple" distributions, see 
for example [32j, [40j or [Ij. For more "complicated" distributions alternative techniques are 
usually employed such as root finding or interpolation. As mentioned above, these methods 
have severe limitations, such as slow execution or setup times. An even bigger setback with 
these techniques however is their failure in the tail regions. The goal of this section is to 
devise an algorithm to overcome these problems. 

Ideally one would like to construct the best rational approximation R*^ „ {v) to the quan- 
tile function Q in the minimax sense. There are numerous methods available such as the 
second algorithm of Renies [36j for the construction R*^ ^ {v). However these methods require 
several a large number of function evaluations, and since Q may only be evaluated to high 
precision through a slow root finding scheme such algorithms usually lead to unacceptable 
setup times. 

Therefore we suggest four alternative algorithms based on certain series acceleration tech- 
niques Q|. The procedures assume the availability of an integration and root finding routine 
to compute the distribution and quantile functions respectively to full machine precision. 
These routines will be used to compute the initial conditions and manage the error in the 
approximation of Q. The inputs to the algorithm are: 1) the distribution parameters, 2) the 
required accuracy e. We partition the unit interval and for convenience name each part as 
follows, 



(0, tl) = Left Tail Region 

[ti,Ui) = Left Region 

[mi,M2] = Central Region 

{u2,Tfi] = Right Region 

(1 - r/j, 1) = Right Tail Region 

The idea behind the first algorithm is then as follows, 

• The asymptotic expansions of Q as developed above are employed to approximate 
Q in the left and right tail regions. Of course these asymptotic expansions are diver- 
gent; however we have found constructing Fade approximants and Levin-type sequence 
transforms particularly useful in summing the series. The order of the approximant is 



^ All of which have been prototyped in Mathematica using double precision arithmetic, and do not rely 
on any of Mathematica's symbolic capabilities, making them portable to lower level languages such as C++ 
or Fortran. 
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chosen in advance, and a numerical search is conducted to determine Tl and r^j, which 
are typically small values f^ 10~^. 

• We will choose the mode quantile location Um '■= F (x^), where Xm is the mode of the 
distribution as the midpoint of the central region. On this region the Taylor expansion 
of Q at Um serves as a useful approximation. Hence we construct the sequence of 
corresponding Pade approximants along the main diagonal of the Pade table. The 
sequence is terminated when a convergent is found which satisfies the required accuracy 
goals. Note like Taylor polynomials, Pade approximants provide exceptionally good 
approximations near a point, in this case the point of expansion p^, but the error 
deteriorates as we move away from the point. Hence we need only check the accuracy 
requirements are met at the end points of the interval [^1,^2] Q We will discuss how 
to choose the points Ui and U2 below. 

• On the left and right regions the left and right solutions of the recycling equation (11. 5p 
denoted A^ and Aft respectively, serve as particularly good approximations. Note that 
this is precisely what they were designed to do. Again a series acceleration technique 
such as Levin's u-transform [38j or Pade summation may be applied to the Taylor 
polynomials of Al and A^^. As is usual with such techniques we observed analytic 
continuation and increased rates of convergence using these techniques. The points at 
which we impose the initial conditions are critical. For now we have chosen uq = Fb {zq) 
where zq = [Qb [tl) + Qb (tim/2)) /2 and zq = [Qb (1 - tr) + Qb [um/'^)) /2 for the 
left and right problems respectively. But by varying these initial conditions we alter the 
range of distribution parameters for which the algorithm is valid. An optimal choice 
has yet to be set. 

The points ui and U2 enclosing the central region are determined by an estimate f of the 
radius of convergence r for the series. 



Q{u) = J2(ln{u-uor. {6A] 

n=0 

In particular we set, 

Ui=Pm- f, 
and, 

U2=Pm + r, 

where, 

r = min{f, \pm - 0.1| , \pm - 0.9|} 



■^ This fact is not entirely true, since in some rare cases a Pade approximant may exhibit spurious poles 
not present in the original function Q. A more robust algorithm would check that none of the real roots of 
the denominator polynomial appearing in the Pade approximant lie in the interval [ui,U2]. Such poles are 
called defects, see [26J for details. 
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For simplicity our approximation of f will be based on the Caucliy-Hadamard formula |22[ 
§ 2.2]. However note that the problem of estimating the radius of convergence is a rather old 
one and many more advanced techniques have been developed to determine f. For example 
Chang and Corliss [9j form a small system of equations based on three, four and five terms 
of the series (16. ip to determine f. However we are not overly concerned if f over estimates 
the radius of convergence of 16. II since this will be compensated by the fact that applying an 
appropriate summation technique will provide analytic continuation of (16. ip . Thus for this 
iteration of the algorithm we will be content with a simple estimate of f provided by the 
Cauchy-Hadamard formula. 

We demonstrate the performance of the algorithm by observing a few test cases for the 
hyperbolic distribution. Note that it would be difficult to test the validity of the algorithm 
for the entire parameter space of the distribution due to the recursive nature of the coeffi- 
cients appearing in the series expansions, so we have biased our testing to parameters which 
frequently occur in financial data and a few extreme cases. As a test case consider the pa- 
rameters a = 89.72, /3 = 4.7184, 5 = 0.0014 and /i = -0.0015. These are the estimated 
parameters for BMW stock prices reported in p2J, computed over a 3 year period. Setting 
the required accuracy to e = 2.98 x 10~®, the algorithm constructs an approximant within 
0.32 seconds. This time is for a rather unoptimized Mathematica prototype of the algorithm 
running on an Intel 17 laptop. We would expect production code to be a fraction of this time. 
The resulting error plots are given in figure 16. 1[ Since the highest degree of the approximant 
is only 24 we expect the algorithm to have reasonably fast execution times. 

Despite meeting the accuracy requirements for the given set of parameters, as can be 
seen from figure 16.11 the error curve produced by this algorithm is far from optimal in the 
minimax sense. Thus our next two algorithms have to do with constructing so called near 
best approximants, which are often good enough in practice due to the difficulties involved 
in finding R^^- To this end consider the partition tl < Ui < U2 < 1 — Tf>. For simplicity we 
choose tl = Tfi = 10~^°. We are interested in the problem of constructing the Chebyshev 
series expansions of the functions Al{z{u)), Q {u) and A^{z{u)) restricted to the sets 
[t,u], [ui,U2] and [m2, 1 — t^] respectively. In each case we introduce a linear change of 
variable x which maps the restricted domain onto the set [—1,1]. To ease notation let 
g (x) := Q|u^ U2] (^)' ^^ ^^^^ '^^^^ ^1 [■^ij'^s] — >■ [— 17 1] is defined by, 

x(.) = !i^ii^^ = !i^^. (6.2) 

2(^2 -Ml) r 

Following from the properties of Q, the function g is continuous and of bounded total vari- 
ation and thus admits the expansion, 

_ 00 
9{x) = ^ + Y,~9i'Tk{x), xe(0,l), (6.3) 

where T^ (x) are the Chebyshev polynomials of the first kind and the coefficients c^ are 
defined by, 

9k = - 9{x) . dx. (6.4) 

vr y_i VI - x^ 
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Since the quantile functions Q we have considered in this report are assumed to be infinitely 
differentiable, elementary Fourier theory tells us that the error made in truncating the series 
f l6.3p after K terms goes to zero more rapidly than any finite power of 1/K as K ^ oo, [TH| 
§ 3]. Such rapid decrease of the remainder motivates us to seek efficient methods to evaluate 
the integral in (16 ■4p . This task can rarely be performed analytically so the the usual process 
is to write (16 ■4p as, 

2 /"'" 
gk = — 9 (cos 6) cos kOdO, 
^ Jo 

and apply a variant of the fast Fourier transform |2L § 29.3]. This method however would 
require several thousand evaluations of the function g and result in unacceptable setup 
times. Fortunately an alternative approach is provided by Thacher |35]. Let the Taylor 
series expansion of g be given by, 

oo 

g{x) = J29kX^ (6.5) 

fe=0 

Then by inverting (16. 2p and substituting into (16. ip we see that the coefficients appearing in 
(16. 5 p may be written as g^ = fq^. Now substituting the following relationship between the 
monomials x^ and the Chebyshev polynomials [ISl eq. 4], 

k 

x''=To{x) + J2ekjT^{x), 

where, 

(2^-U,i,) j-keven 
[O i-koM 

into (16. 5p . one may express the Chebyshev coefficients g^ in terms of the Taylor coefficients 
(yffc as follows, 

oo 

9k = '^gjdj,k, (6.6) 

j=k 

Thacher's approach was simple, he observed that applying Shank's transform to (16. 6p one 
may approximate the Chebyshev coefficients even for slowly convergent series. In our ex- 
perimentation we found Levin's w-transform to also be affective, which of course was not 
discovered at the time Thacher wrote his paper. Note that Levin's u transform has been 
reported in many instances to outperform Shank's transform [31]. For efficient Fortran and 
C++ implementations of these two transforms see |^46j and [35j respectively. Thus given the 
knowledge of the Taylor coefficients g^. we now have a method to approximate the Chebyshev 
coefficients without having to resort to the computationally expensive task of evaluating (16. 4p 
through numerical quadrature. An estimate of the error made in truncating the series (16.30 
is given by the magnitude of the first neglected term. In a similar fashion we can construct 
Chebyshev series expansions for the restrictions Ai|r ^ , and ^hIu^ i-r i- ^^ ^^^ ^^ ^^^^ 



6 Numerical Techniques and Examples 



from figure 16.21 an algorithm based on constructing truncated Chebyshev expansions in this 
manner provides us with a much more satisfactory error curve. 

Perhaps the greatest drawback of this algorithm is that it requires one to compute a large 
number of Chebyshev coefficients. In our example the series needed to be truncated at the 
35**^ term in order to satisfy the accuracy requirements in the right region, see figure [6721 We 
stress that approximating these 35 coefficients via l6.6l is not a computationally expensive task 
in contrast to a numerical integration scheme, however evaluating a 35**^ degree polynomial 
may lead to unacceptable execution times. It is for this reason we consider constructing 
Chebyshev- Pade approximants of the function g, which are rational functions of the form, 

rj. ( \ ._ Qo?o (a;) H h amTm jx) 

boTo (x) H h bnTn (x) 

which have a formal Chebyshev series expansion in agreement with (16. 3p up to and including 
the term Qm+nTm+n (x). Given the knowledge of the first n + m + 1 Chebyshev coefficients 
^0, • • • iQn+mi an efficient way to construct Tm,n {x) is to employ the algorithm of Sidi ^, 
in which the coefficients appearing in 16.71 are computed recursively. For an estimate of the 
error given by (16.71) see [371, ^q. 7.6-12]. Applying this scheme to the same test case as above 
we obtain the error plots given in figure ( 16. 3p . 

The real benefit of this algorithm is the fast setup and execution times it provides. For 
example we required a truncated Chebyshev series of degree 34 to approximate the hyperbolic 
quantile to the required accuracy in the left region of our test case, again see figure ( 16. 2p . 
On the other hand a Chebyshev-Pade approximant of degree only 8 provides us with a much 
better approximant in the same region, contrast with figure (16. 3p . 

Our final algorithm is based on multipoint Pade approximants [26]. Consider the partition 
Tl < Pm < 1 — tr, where Pm is the mode location. Choose the points Uq, . . . , Uk G [TL,Pm\, 
and let So, . . . , s/^- G N. Suppose now we want to build a rational approximant Rm,n valid on 
[TLiPm\ which satisfies the following conditions, 

Af [zk) = R'ilAzk) , 0<k<K,0<j<Sk-l, (6.8) 

where Zk = Qb (uk)- Such a problem is called an osculatory rational interpolation problem. 
The interpolation data on the right hand side of ( 16. 8 P may be generated by solving the 
recycling equation (11. 5p for A^ with initial conditions imposed at the points Uq, . . . , Uk- 
Provided a solution to the problem exists it may be solved efficiently through the generalized 
Q.D. algorithm [31]. This algorithm generates the partial numerators and denominators 
of a continued fraction, who's convergents are rational functions called multipoint Pade 
approximants satisfying (16. 8p . Similarly by considering solutions to A^ we may find a rational 
approximant valid on [pm, 1 — t/j]. By choosing only 3 interpolation points, we apply the 
method to our test case and provide the results in figure (16. 4|) . 

One could go a step further and attempt to construct the minimax approximant -R^n) 
however as mentioned most of the methods known to the present author to construct R*^ ^ 
require evaluating Q at a large number of points. It is for this reason we recommend Maehly's 
indirect method [31j . This method of constructing R^^^ i^ applicable when a sufficiently 
accurate approximation of Q is available. Such an initial approximation may be obtained by 
the methods described above. Maehly's indirect method then seeks to find a more optimal 
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approximation in the Chebyshev sense without the need to evaluate Q, which of course in 
our case is an expensive operation. Thus one would expect reduced setup up times with this 
approach, however we have not explored this idea further. 
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